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Analysis  of  Different  Approaches  to  Modeling  of  Nozzle 
Flows  in  the  Near  Continuum  Regime  (Preprint) 

E.  V.  Titov,  Rakesh  Kumar,  D.  A.  Levin*  and  N.  E.  Gimelshein,  S.  F.  Gimelshein' 
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Abstract.  Conical  nozzle  flows  are  studied  for  Reynolds  numbers  of  1,230  and  12,300  using  different  numerical  techniques: 
DSMC  Method,  Navier-Stokes/CFD  accounting  for  velocity  slip  and  temperature  jump  boundary  conditions,  and  statistical 
and  deterministic  approaches  to  the  solution  of  BGK  equation.  Detailed  comparisons  of  the  stability,  accuracy,  and  conver¬ 
gence  of  the  employed  numerical  techniques  provides  better  understanding  of  their  benefits  and  deficiencies,  and  assists  in 
selecting  the  most  appropriate  technique  for  a  particular  nozzle  and  flow  application.  The  deterministic  and  statistical  solutions 
of  the  BGK  equation  were  found  to  be  in  good  agreement  with  the  benchmark  DSMC  results.  The  Navier-Stokes  solution 
differs  from  DSMC  in  the  boundary  layer. 
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INTRODUCTION 

Accurate  and  numerically  efficient  modeling  of  low 
and  moderate  Reynolds  number  nozzle  flows  is  of¬ 
ten  problematic  due  to  the  existence  of  multiple  flow 
length  scales.  Experimental  studies  of  micro-nozzle 
based  thrusters  are  rare,  expensive,  and  may  not  provide 
the  necessary  precision  in  the  measurement  of  principal 
nozzle  characteristics  such  as  thrust,  flow  rate,  and  spe¬ 
cific  impulse.  The  problem  becomes  even  more  severe 
when  the  influence  of  MEMS  based  thrusters  on  sensi¬ 
tive  spacecraft  surfaces  needs  to  be  analyzed.  Back  flow 
produced  by  such  devices  plays  a  major  role  in  the  con¬ 
tamination  of  sensitive  electronic  devices  such  as  optical 
instruments,  and  solar  panels  etc.,  which  in  turn  may  ad¬ 
versely  affect  the  life  span  of  spacecraft.  This  back  flow 
formation  is  sensitive  to  the  flow  conditions  at  the  nozzle 
lip,  and  is  also  known  to  be  difficult  to  study  experimen¬ 
tally  for  such  micro-nozzle  flows  [1]. 

The  development  of  accurate  numerical  tools  capable 
of  handling  micronozzle  flows  is  therefore  important,  but 
challenging  at  the  same  time  because  the  flow  regime 
changes  from  continuum,  near  the  nozzle  throat,  to  tran¬ 
sitional  at  the  nozzle  exit.  Both  kinetic  methods,  such  as 
the  Direct  Simulation  Monte  Carlo  (DSMC),  and  con¬ 
tinuum,  techniques  based  on  Navier-Stokes  (N-S)  solu¬ 
tions,  must  meet  computational  and  physical  challenges 
when  applied  to  these  flows.  The  major  problem  with  the 
DSMC  method  [2]  is  the  associated  computational  cost 
when  high  density  portions  of  the  flow  have  to  be  accu¬ 
rately  modeled.  On  the  other  hand,  conventional  contin¬ 
uum  CFD  techniques  are  inapplicable  in  the  regions  of 
high  gradients  and  strong  rarefaction  even  with  the  use 


of  velocity  slip  and  temperature  jump  boundary  condi¬ 
tions  at  the  nozzle  surface.  It  would  therefore  be  highly 
desirable  to  have  a  single  method  that  allows  an  accurate 
and  efficient  one- step  modeling  of  high  density  nozzle 
and  low  density  plume  flows.  To  this  end,  the  use  of  sim¬ 
plified  forms  of  the  Boltzmann  equation,  usually  called 
model  kinetic  equations,  such  as  Bhatnagar-Gross-Krook 
(BGK)  may  be  useful  [3] .  A  method  based  on  BGK  equa¬ 
tion  is  expected  to  be  more  efficient  than  the  DSMC 
method  in  the  continuum  and  near-continuum  regime, 
and  more  accurate  than  the  solution  of  the  N-S  equations 
in  the  transition  regime. 

Recent  years  have  witnessed  a  renewed  interest  and 
significant  advances  in  the  solution  of  model  kinetic 
equations  such  as  BGK,  with  deterministic,  either  finite 
difference  or  finite  volume,  approaches  typically  used 
in  the  solution  procedure.  A  particle  approach  to  obtain 
the  solution  to  the  BGK  equation  was  first  proposed  in 
Ref.  [4].  It  was  then  extended  to  model  the  ES-BGK 
equation  in  Ref.  [5],  and  further  extended  to  include  ro¬ 
tational  degrees  of  freedom  in  Ref.  [6].  The  main  goal 
of  this  paper  is  to  apply  statistical  and  deterministic  ap¬ 
proaches  to  obtain  the  solution  of  the  BGK  equation  to 
simulate  rarefied  gas  flow  through  a  conical  nozzle,  and 
compare  the  results  in  terms  of  accuracy  and  efficiency  to 
the  solutions  obtained  with  the  traditional  DSMC  method 
and  the  Navier-Stokes  equations.  Four  different  numeri¬ 
cal  approaches  -  DSMC,  finite  volume  and  particle  solu¬ 
tions  of  the  BGK  and  ES-BGK  model  kinetic  equations 
and  an  equilibrium  DSMC  (eDSMC)  technique  are  used 
to  study  nozzle  flows  expanding  into  a  vacuum. 
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GEOMETRY  AND  FLOW  CONDITIONS 

Gas  flow  of  pure  argon  through  a  conical  nozzle  into 
vacuum  is  considered  in  this  work.  The  diverging  part 
of  the  nozzle  is  modeled,  and  its  geometry  is  taken  from 
Ref.  [7].  The  nozzle  throat  diameter  is  2.5  mm,  the  length 
of  the  diverging  part  is  50.7  mm,  and  the  half-angle  is 
20deg.  The  surface  temperature  of  the  nozzle  is  assumed 
to  be  300  K.  Numerical  results  are  obtained  for  two 
throat-diameter  based  Reynolds  numbers,  1,230  (Case  I) 
and  12,300  (Case  II),  with  the  stagnation  temperature  of 
333  K.  In  all  numerical  approaches,  the  computational 
domain  starts  at  the  nozzle  throat  and  covers  the  entire 
diverging  part  of  the  nozzle,  as  well  as  a  small  part  of  the 
plume  to  avoid  the  influence  of  the  downstream  boundary 
conditions  [8]. 

COMPUTATIONAL  METHODS 
DSMC  method 

The  DSMC-based  SMILE  computational  solver  was 
used  in  this  work.  Details  on  the  tool  may  be  found 
in  Ref.  [9].  The  SMILE  features  that  were  used  in  the 
present  work  include  axisymmetric  capability  with  radial 
weights,  different  grids  for  collisions  and  macroparam¬ 
eters,  both  of  which  are  two-level  adaptable  Cartesian 
grids,  and  parallel  implementation  with  efficient  load 
balancing  techniques.  The  majorant  frequency  scheme 
was  employed  for  modeling  molecular  collisions  in  use 
of  DSMC  [10].  The  VHS  model  was  used  for  modeling 
molecular  interactions.  Diffuse  reflection  with  full  ther¬ 
mal  accommodation  was  assumed  on  the  nozzle  wall. 
For  both  Reynolds  numbers,  solutions  independent  of 
grid,  time  step,  and  number  of  particles  were  obtained. 
For  Re=  1,230,  there  was  virtually  no  difference  observed 
between  solutions  obtained  for  0.3  million  cells  with  1 
million  molecules,  and  3  million  cells  with  10  million 
molecules.  For  Re=  12,300,  this  was  true  for  numerical 
parameters  up  to  an  order  of  magnitude  larger.  Tables  1 
and  2  provide  a  summary  of  the  numerical  parameters  for 
the  two  cases  considered  in  this  work. 


Solution  of  NS  equations 

A  commercial  code,  CFD++  [11]  has  been  used  in 
this  work  to  solve  the  Navier- Stokes  equations.  CFD++ 
is  a  flexible  CFD  software  suitable  for  the  solution 
of  steady/unsteady,  compressible/incompressible  N-S 
equations,  including  multi- species  capability  for  perfect 
and  reacting  gases.  In  this  work,  a  perfect-gas  compress¬ 
ible  N-S  solver  was  used  with  second  order  spatial  dis¬ 


cretization  and  implicit  time  integration.  Second  order 
velocity  slip  and  temperature  jump  conditions  were  im¬ 
posed  on  the  nozzle  wall.  A  supersonic  inflow  with  pre¬ 
scribed  parameters  was  applied  at  the  nozzle  throat,  and 
backpressure  of  1  Pa  was  imposed  at  the  outflow  bound¬ 
aries.  The  results  presented  in  this  paper  were  obtained 
for  a  multi-block  rectangular  grid  with  a  total  of  14,400 
nodes.  The  computations  were  also  conducted  for  four 
times  smaller  and  four  times  larger  numbers  of  nodes, 
and  found  to  be  fully  grid-resolved  with  14,400  nodes. 

Finite  Volume  method  for  BGK  and 
ES-BGK  equations 

A  finite  volume  solver  SMOKE  developed  at  ERC  has 
been  used  to  deterministically  solve  the  BGK  and  ES- 
BGK  equations.  SMOKE  is  a  parallel  code  based  on  con¬ 
servative  numerical  schemes  developed  by  L.  Mieussens 
[12].  A  second  order  spatial  discretization  with  axial 
symmetry  is  used  along  with  implicit  time  integration. 
A  supersonic  inflow  condition  is  used  at  the  nozzle 
throat,  and  vacuum  outflow  conditions  are  set  at  the  outer 
boundaries.  Fully  diffuse  reflection  with  complete  en¬ 
ergy  accommodation  is  applied  at  the  nozzle  surface. 
The  spatial  grid  convergence  was  achieved  by  increas¬ 
ing  the  number  of  nodes  from  3,600  to  over  17,000. 
The  convergence  on  the  velocity  grid  is  also  studied,  with 
the  number  of  (x,  r,  0)  points  ranging  from  (20,10,18)  to 
(30,35,50). 


Statistical  method  for  BGK  equation 

In  our  earlier  work  [13],  we  had  developed  and  stud¬ 
ied  a  statistical  technique,  called  eDSMC,  which  models 
continuum  flows  through  a  collision  enforcement  proce¬ 
dure  which  guarantees  full  relaxation  of  the  molecular 
thermal  velocities  to  the  state  of  local  equilibrium.  The 
technique  is  able  to  solve  inviscid  flows  with  tangency 
boundary  conditions  at  the  wall,  but  tends  to  under  pre¬ 
dict  the  viscous  effects  in  the  boundary  layer  when  dif- 
fusional  boundary  conditions  are  used.  In  this  work  we 
continue  our  effort  to  apply  statistical  methods  to  mod¬ 
erate  and  high  Reynolds  number  nozzle  flows  by  making 
use  of  the  BGK  [3]  model. 

Recently,  a  number  of  authors  [5,  14,  6]  have  devel¬ 
oped  particle  approaches  to  the  solution  of  the  BGK  and 
ES-BGK  model  equations.  The  essence  of  these  kinetic 
approches  is  to  relax  the  flow  to  local  Maxwellian  or 
Gaussian  equilibrium  by  choosing  a  fraction  of  simulated 
particles  available  in  a  computational  cell,  and  assigning 
them  new  velocities  according  to  the  local  Maxwellian 
(or  ellipsoidal)  distribution.  If  the  collision  frequency  for 
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such  a  velocity  reassignment  is  properly  computed  and 
local  values  of  the  translational  temperature  in  the  cell 
are  known,  then  the  procedure  should  simulate  the  colli¬ 
sion  term  on  the  right  hand  side  of  the  BGK  equation, 

d  ( nf )  jdt  +  c-d  ( nf )  / dr  =vn(fe-f ) 

where  n  is  the  number  density,  V  is  the  collision  fre¬ 
quency,  /  is  the  molecular  velocity  distribution  function, 
fe  is  the  Maxwellian  distribution  function,  c  is  the  veloc¬ 
ity  vector,  and  r  is  the  position  vector. 

The  details  of  the  statistical  BGK  scheme  are  as  fol¬ 
lows.  The  collision  frequency  is  calculated  as 

/  T®  \ 

V  =  Pr-nk  (  — — )  Tl~m 

where  Pr  is  the  Prandtl  number  (1  for  the  BGK  equa¬ 
tion),  k  is  the  Boltzmann  constant,  T  is  the  translational 
temperature  in  the  flow,  (iref  is  the  gas  dynamic  viscos¬ 
ity  at  Tref.  The  collision  frequency  is  calculated  for  each 
computational  cell  at  each  time  step  based  on  the  local 
translational  temperature  T  and  the  local  number  density 
n.  The  local  number  density  n  was  averaged  over  a  large 
number  of  computational  time  steps,  while  the  local  tem¬ 
perature  T  is  computed  based  on  the  instantaneous  ther¬ 
mal  velocities  of  the  computational  particles  in  the  cell. 

The  number  of  particles  in  a  cell  preselected  for  veloc¬ 
ity  re- sampling  was  calculated  as  follows: 

Nc  =  int(V(l  —  exp(— vAf)) 

where  At  is  the  time  step,  N  is  the  number  of  particles 
in  the  cell  and  int  operator  means  the  nearest  smaller 
integer.  In  order  to  compensate  for  the  systematic  error 
that  such  an  operator  produces,  one  more  particle  was 
added  to  the  list  of  preselected  particles  with  the  proba¬ 
bility  given  by  following  equation: 

PC=N(  1  —  exp(— vAt)  —  int(V(l  —  exp(— vAt)) 

The  preselected  particles  receive  new  velocities  ac¬ 
cording  to  a  Maxwellian  distribution  corresponding  to 
the  local  cell  temperature  and  velocity.  The  velocities 
of  particles  which  have  not  been  preselected  remain  un¬ 
changed  in  the  current  timestep. 

Although  the  technique  is  not  limited  to  simple 
gases  [6],  argon  was  used  as  the  working  gas  in  or¬ 
der  to  understand  the  basic  features  and  limitations  of 
the  method,  avoiding  the  additional  complication  of 
translational-rotational  relaxation.  In  order  to  reduce  the 
statistical  error  so  that  small  differences  among  the  dif¬ 
ferent  methods  may  be  studied,  a  large  number  of  simu¬ 
lated  particles  were  used  such  that  the  minimum  number 
of  particles  per  cell  was  about  150.  This  unusually  large 
number  of  particles  per  cell  also  provided  a  better  sta¬ 
tistical  approximation  of  local  instantaneous  cell  based 
temperature. 


RESULTS  AND  DISCUSSION 
Low  pressure  Case  I 

The  results  presented  below  for  the  low  pressure  Case 
I  were  obtained  by  the  DSMC  code  (SMILE),  its  mod¬ 
ification  implementing  the  statistical  BGK  scheme,  a  fi¬ 
nite  volume  (FV)  ES-BGK  solver,  and  a  Navier- Stokes 
solver,  CFD++.  To  quantify  the  differences  among  the 
different  approaches,  the  gas  parameters  are  presented 
along  two  cross  sections,  the  nozzle  centerline  and  the 
nozzle  exit  plane.  The  velocity,  and  temperature  along 
the  centerline  are  shown  in  Figs.  1  and  2.  Here,  X=0  cor¬ 
responds  to  the  nozzle  throat.  It  is  evident  that  all  four 
of  our  approaches  produce  very  similar  results  for  all  gas 
properties  under  consideration. 

The  axial  velocity  and  temperature  along  the  nozzle 
exit  plane  are  presented  in  Figs.  3  and  4.  Here,  Y=0  cor¬ 
responds  to  the  nozzle  centerline.  The  FV  ES-BGK  so¬ 
lution  is  close  to  that  of  DSMC  for  both  the  axial  veloc¬ 
ity  and  gas  temperature.  The  particle  BGK  result  agrees 
well  with  the  FV  ES-BGK  and  DSMC.  The  NS  solver 
produces  results  accurate  in  the  coreflow  but  strongly 
different  from  the  ES-BGK,  DSMC  and  particle  BGK 
throughout  the  boundary  layer.  Generally,  we  can  con¬ 
clude  that  the  statistical  and  FV  ES-BGK  results  are  in 
very  good  agreement  with  the  DSMC  in  the  entire  diverg¬ 
ing  part  of  the  nozzle.  The  solution  of  the  model  kinetic 
equation  accurately  captures  both  the  boundary  layer  and 
the  nozzle  coreflow. 


High  pressure  Case  II 

Results  for  the  high  pressure  Case  II  were  obtained  by 
the  baseline  DSMC,  statistical  BGK  scheme,  Finite  Vol¬ 
ume  BGK  solver,  the  Navier- Stokesi  solver  and  eDSMC 
method.  A  quantitative  comparison  of  flow  velocity  and 
temperature  for  Case  II  is  given  in  Figures.  5  and  6, 
where  the  velocity  and  temperature  profiles  along  the 
nozzle  centerline  are  shown,  and  in  Figures.  7  and  8  for 
the  velocity  and  temperature  profiles  along  the  nozzle 
exit  plane.  It  is  clearly  seen  that  the  DSMC,  NS,  FV  ES- 
BGK  and  statistical  BGK  profiles  of  both  velocity  and 
temperature  are  very  close  along  the  nozzle  centerline. 
The  results  obtained  with  eDSMC  are  similar  to  parti¬ 
cle  BGK  in  the  coreflow.  The  results  along  the  exit  plane 
show  that  there  is  still  some  difference  between  the  NS 
and  DSMC  predictions,  both  in  temperature  and  flow  ve¬ 
locity.  As  in  Case  I,  this  difference  is  explained  by  the 
limitations  of  the  velocity  slip  boundary  condition  used 
in  NS.  There  is  some  impact  of  flow  nonequilibrium  too, 
but  its  contribution  is  significant  only  near  the  nozzle  lip, 
where  Tx/T  ratio  reaches  0.85. 
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Tables  1  and  2  show  the  comparison  of  computational 
time  required  by  various  methods  for  low  pressure  and 
high  pressure  cases  respectively.  It  is  interesting  to  note 
that  computational  efforts  for  the  DSMC  and  statistical 
BGK  schemes  were  of  the  same  order,  but  the  conver¬ 
gence  process  was  quite  different.  The  DSMC  schemes 
tend  to  converge  faster,  but  requires  more  computational 
effort  to  collect  the  sufficient  information  for  the  solution 
to  be  smooth.  On  the  other  hand  statistical  BGK  con¬ 
verges  slower  (perhaps  due  to  the  "history"  of  the  macro- 
porameter  sampling  procedure  which  defines  the  local 
Maxwellian  distribution  function),  but  the  smoothness  of 
the  results  is  achieved  earlier  than  in  the  case  of  DSMC 
(probably  for  the  same  reason).  eDSMC,  although  very 
fast  (comparable  with  the  NS  solver),  does  not  capture 
the  viscous  effects  in  the  boundary  layer,  yet  solution  in 
the  inviscid  core  of  the  flow  is  remarkably  close  to  the 
solutions  obtained  by  using  other  methods. 

CONCLUSIONS 

Argon  flow  through  a  conical  nozzle  was  studied  for 
two  Reynolds  numbers  of  1,230  and  12,300,  using  four 
different  approaches.  These  include  one  continuum  ap¬ 
proach  (solution  of  Navier-Stokes  equations)  and  three 
kinetic  approaches  (the  DSMC  method,  and  statistical 
and  deterministic  methods  for  the  BGK  equation).  Anal¬ 
yses  of  the  accuracy  of  the  approaches  and  their  numer¬ 
ical  efficiency  was  conducted.  Several  conclusions  can 
be  drawn  from  the  results  of  the  computations.  The  most 
accurate  data  in  both  high  and  low  pressure  cases  was  ob¬ 
tained  by  DSMC  even  though  parameters  of  the  numer¬ 
ical  scheme  were  somewhat  relaxed  in  the  high  pressure 
Case  II. 

The  Navier-Stokes  solutions  are  in  good  agreement 
with  the  DSMC  results  in  the  higher  density  portion  of 
the  flow  and  in  the  coreflow,  where  rarefaction  effects  are 
small.  In  the  boundary  layer,  even  though  the  velocity 
slip  and  temperature  jump  boundary  conditions  were 
used,  there  is  a  difference  between  the  NS  and  DSMC 
solutions.  The  statistical  and  finite  volume  solution  of 
the  ES-BGK  equation  are  in  fair  agreement  with  the 
DSMC  method  in  the  entire  computational  domain  for 
both  Reynolds  numbers. 
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FIGURE  1.  Velocity  profile  along  nozzle  axis  for  low  pres¬ 
sure  case 


FIGURE  3.  Velocity  profile  at  nozzle  exit  plane  for  low 
pressure  case 


FIGURE  2.  Temperature  profile  along  nozzle  axis  for  low 
pressure  case 


FIGURE  4.  Temperature  profile  at  nozzle  exit  plane  for  low 
pressure  case 
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TABLE  1.  Parameters  of  the  methods,  Case  I 


Method 

DSMC 

Part  BGK 

FV  ES-BGK 

NS 

Time  (CPUH) 

200 

1000 

400 

<1 

Number  of  particles 

10  mil 

50  mil 

- 

- 

Number  of  cells 

3  mil 

0.5  mil 

3,600 

3,600 

TABLE  2.  Parameters  of  the  methods,  Case  II 


Method 

DSMC 

Part  BGK 

FV  ES-BGK 

NS 

eDSMC 

Time  (CPUH) 

1000 

1000 

5000 

3 

100 

Number  of  particles 

50  mil 

50  mil 

- 

- 

1  mil 

Number  of  cells 

15  mil 

0.5  mil 

17,200 

14,400 

100,000 

FIGURE  5.  Velocity  profile  along  nozzle  axis  for  high  pres-  FIGURE  7.  Velocity  profile  at  nozzle  exit  plane  for  high 
sure  case  pressure  case 


FIGURE  6.  Temperature  profile  along  nozzle  axis  for  high  FIGURE  8.  Temperature  profile  at  nozzle  exit  plane  for  high 
pressure  case  pressure  case 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


